Cytoskeletal Keratins Are Overexpressed in a Zebrafish Model of Idiopathic Scoliosis

Idiopathic scoliosis (IS) is a three-dimensional rotation of the spine >10 degrees with an unknown etiology. Our laboratory established a late-onset IS model in zebrafish (Danio rerio) containing a deletion in kif7. A total of 25% of kif7co63/co63 zebrafish develop spinal curvatures and are otherwise developmentally normal, although the molecular mechanisms underlying the scoliosis are unknown. To define transcripts associated with scoliosis in this model, we performed bulk mRNA sequencing on 6 weeks past fertilization (wpf) kif7co63/co63 zebrafish with and without scoliosis. Additionally, we sequenced kif7co63/co63, kif7co63/+, and AB zebrafish (n = 3 per genotype). Sequencing reads were aligned to the GRCz11 genome and FPKM values were calculated. Differences between groups were calculated for each transcript by the t-test. Principal component analysis showed that transcriptomes clustered by sample age and genotype. kif7 mRNA was mildly reduced in both homozygous and heterozygous zebrafish compared to AB. Sonic hedgehog target genes were upregulated in kif7co63/co63 zebrafish over AB, but no difference was detected between scoliotic and non-scoliotic mutants. The top upregulated genes in scoliotic zebrafish were cytoskeletal keratins. Pankeratin staining of 6 wpf scoliotic and non-scoliotic kif7co63/co63 zebrafish showed increased keratin levels within the zebrafish musculature and intervertebral disc (IVD). Keratins are major components of the embryonic notochord, and aberrant keratin expression has been associated with intervertebral disc degeneration (IVDD) in both zebrafish and humans. The role of increased keratin accumulation as a molecular mechanism associated with the onset of scoliosis warrants further study.

Over the last decade, D. rerio (zebrafish) has emerged as a popular model for human scoliosis based on similar vertebral anatomy, ease of breeding, and genetic manipulation [21]. To date, multiple zebrafish genetic mutants (e.g., ptk7 [22,23], POC5 [24], kif6 [25]) have linked genes involved in cilia and ciliary motility to a progressive scoliosis, giving potential insight into molecular mechanisms that may regulate scoliosis development.
These mutants exhibit ciliary motility defects and often develop hydrocephalus and/or a disruption of the Reissner fiber, a fiber aggregation along the central canal that is essential to zebrafish axial development [26][27][28]. Mechanisms related to cerebrospinal fluid (CSF) flow [23,29], specialized contacting neurons within the axial sensory system (CSF-cNs) [30], inflammation [31], and urotensin peptides [32] have been hypothesized as related to these phenotypic observations; however, a definitive mechanism related specifically to scoliosis development has not been characterized.
Our zebrafish model for scoliosis, kif7 co63/co63 , develops juvenile-onset spinal curvatures without obvious vertebral malformations, hydrocephalus, or malformations in the Reissner fiber [33]. As patients with IS are not yet known to have gross morphological changes to these structures and develop scoliosis during the juvenile-adolescent period, kif7 co63/co63 zebrafish may uniquely mirror the human condition. KIF7 encodes a broadly conserved kinesin ciliary protein that localizes to the axonemal tip of primary cilia and binds to the plus-ends of microtubules [34]. The protein functions as a scaffold protein for ciliary function and acts as both a negative and positive regulator of the hedgehog (Hh) signaling pathway, an evolutionary conserved molecular pathway central to embryonic development, limb patterning and musculoskeletal maintenance [34], including formation and maintenance of the intervertebral disc [35]. Kif7 primarily acts by suppressing the Gli1 transcription factor. In zebrafish, Kif7 accumulates at the ciliary tip, as observed in mammals, as well as within cytoplasmic puncta, which sequester Gli1 and Gli2 and disperse in response to Hh pathway activation [36]. In humans, loss of function mutations in KIF7 have been linked to both Joubert and the rare acrocallosal syndromes (OMIM #611254), two ciliopathies with overlapping, system-wide defects including developmental disability, skeletal abnormalities and kidney disease. Scoliosis has been observed in 5-33% of Joubert syndrome cases [37][38][39], potentially related to early hypotonia, whereas the scoliosis prevalence in acrocallosal syndrome has not yet been described.
Approximately 25% of kif7 co63/co63 zebrafish develop spinal curvatures as juveniles with no evidence of abnormalities in brain morphology or hydrocephaly, and no morphological changes to the central canal cilia or the Reissner fiber. Our hypothesis is that study of the kif7 co63/co63 zebrafish, both with and without the scoliosis phenotype, will provide insight into potential molecular mechanisms underlying scoliosis development. To this end we performed bulk transcriptome mRNA sequencing analyses of our mutant kif7 co63/co63 embryonic (4 days post fertilization) and young adult (6 weeks past fertilization [wpf]) zebrafish with and without scoliosis, as well as age-matched heterozygous and wild-type controls. Specific findings are validated via qRT-PCR and histology.

Zebrafish Husbandry and Crossing
Zebrafish were housed at the University of Colorado Anschutz Medical Campus aquatic facility. Husbandry and experimental protocols were approved under the Institutional Animal Care and Use Committee Protocol #00370. Animals were maintained in a 14 h light/10 h dark cycle at 28.5 • C. Animals under anesthesia in 168 mg/L Tricaine were fin clipped for genotyping and RNA extraction. Animals were euthanized in 400 mg/L Tricaine, followed by immersion in ice water.
Zebrafish were bred in a sloped 1.7 L static tank. Males and females were placed in the tank separated by a divider. Each tank was set with a maximum of 8 total fish in an either 1:1 or 2:1 female to male ratio. Fish were placed in static tank in the afternoon after their last daily feeding. The following morning, water in the static tanks was changed and the dividers were pulled. Fish were monitored and embryos were harvested by siphoning after two hours.

Genotyping
At 24 h past fertilization (hpf), 8 embryos were euthanized by ice water immersion. Embryos were placed in 50 µL of lysis solution (500 µL 1 M Tris pH 8.3, 2.5 mL 1 M KCL, 1.5 mL 10% Tween, 1.5 mL 10% NP40, 44 mL ddH20). The 6 wpf zebrafish were genotyped by fin clip immediately post euthanasia. Briefly, zebrafish were immersed in 400 mg/L tricaine until gill movement stopped and tail fins were clipped with a sterile razor. Clipped fins were immersed in 50 µL lysis solution and run with the following reaction: 20 min at 95 • C, 2 min on ice, added 2.5 µL Proteinase K (Invitrogen, Waltham, MA, USA), 3 h at 55 • C, 10 min at 95 • C, and hold at 12 • C. Samples underwent the following PCR reaction using 12.5 µL 2× GOTaq Green Master Mix (Promega, Madison, WI, USA), 1 µL gDNA, 0.5 µL 10 µM F and R primers, and 10 µL sterile water (see Table 1): Routine genotyping of zebrafish lines was performed by PCR as above, followed by gel electrophoresis using a 3% agarose gel with a 1:1 ratio of standard agarose to MetaPhor agarose (Lonza, Basel, Switzerland). The gel was run at 120 V for 1 h and 45 min followed by UV imaging (Supplemental Figure S5).

RNA Extraction and Sequencing (Novogene)
At 4 days past fertilization (dpf), 20 embryos per tube were euthanized and pooled in 1 mL of RNA lysis buffer (Zymo Research, Irvine, CA, USA). Embryos were lysed using a VWR Bead Mill Homogenizer with Lysing Matrix A beads (MP Biomedicals, Santa Ana, CA, USA). Samples were iced for 5 min, lysed for 45 s, then iced again for 5 min. This process was completed three times to ensure maximal lysis. The same protocol was used to extract RNA for 6 wpf zebrafish, but 1 fish was used per tube. RNA was extracted using the QuickRNA Miniprep (Zymo) using the on-column DNase digestion protocol.

Bioinformatic Filtering
Bioinformatic filtering was performed by the Novogene Company (Beijing, China) following an established pipeline. Briefly, Illumina image data were base called using CASAVA and stored in FASTQ format. The sequencing error rate is represented in Qphred scores. Low-quality reads (reads where uncertain nucleotides constitute >10% of the total read, or when low-quality nucleotides [base quality < 5] constituted >50% of the read) and reads containing adapters were filtered out. Sequences were mapped to the zebrafish reference genome GRCz11/danRer11 using HISAT2 v.2.0.5. Gene expression levels and differential gene expression analysis were calculated through DESeq2 [41]. Multiple testing correction was conducted using Benjamini and Hochberg's approach for controlling the False Discovery Rate (FDR).
RT-qPCR: Complementary DNA (cDNA) was created from RNA samples using the High-Capacity cDNA Reverse Transcription Kit with RNase inhibitor (Applied Biosystems, Waltham, MA, USA) using random primers. RT-qPCR was conducted using the same RNA as that submitted for RNA sequencing (see RNA Extraction). One microgram of RNA was used per reaction. Real-time quantitative PCR primer sequences are provided in Supplemental File S1 Table S5. Primers were analyzed using IDT's OligoAnalyzer tool (Integrated DNA Technologies, Coralville, IA, USA) to ensure length, GC content, melting temp, hairpin tm, self-dimer and heterodimer parameters are well within range. All primers also checked using BLAT (UCSC Genome Browser, University of California Santa Cruz, Santa Cruz, CA, USA) to ensure specificity. Primer efficiency was tested for all primers to ensure genes of interest and housekeeping gene is comparable. SsoAdvanced SYBR Green Supermix (Bio-Rad, Hercules, CA, USA) was used for reactions according to the manufacturer's instructions. An amount of 2 ng of cDNA was used per reaction. Reactions were run on a Bio-Rad CFX 96 instrument at 95 • C for 10 min followed by 40 cycles of 95 • C for 15 s and 60 • C for 1 min. Results were analyzed in accordance with the 2 −ddCt method [42]. Embryo samples were pooled (n = 20 embryos per pool) and adult samples were run as individuals. We conducted n = 3 biological replicates per run, with all samples in triplicate. All experiments were repeated at least twice by independent technicians.

Keratin Staining and Histology
Tissue staining was performed in accordance to the protocol outlined in [43]. The 6 wpf scoliotic and non-scoliotic fish were fixed with 10% neutral buffered formalin for 3 days at 4 • C. Fish were decalcified using 20% EDTA pH 8 for 10 days at room temperature using nutator. Fish were processed and embedded in paraffin using Tissue-Tek VIP 6 AI and Tissue-Tek TEC. The 6.5 µm sections were cut and dewaxed prior to staining with Mayer's Hematoxylin, Phloxine B, Alcian Blue, and Orange G (Electron Microscopy Sciences, Hatfield, PA, USA, 26401-04, 26401-01, 26401-02, 26401-03). Slides were visualized using Zeiss Axio Scan. Z1 at 20× using brightfield.
Transcriptomes cluster by sample age and genotype: To first determine overall transcriptomic changes occurring within scoliotic versus (vs.) non-scoliotic kif7 co63/co63 zebrafish, we performed bulk RNA sequencing of kif7 co63/co63 , kif7 co63/+ and wild-type (AB) zebrafish at 6 wpf, which is near the time of curvature onset (n = 3 individuals per phenotype). Additionally, we sequenced homozygous, heterozygous, and AB embryos (4 days past fertilization [dpf]; n = 3 embryo pools [20 embryos per pool] per genotype). Sequencing information is provided in Supplemental Table S1. The 6 wpf timepoint was selected for sequencing to correspond to just after the onset of spinal curvature, when scoliosis status is able to be distinguished between individuals [33]. The 4 dpf timepoint was selected due to the expression profile of kif7 (Supplemental Figure S1), which shows an increase in expression at day 5 in kif7 co63/co63 embryos. At 4 dpf, embryos show a reduced expression of kif7 (0.545 fold change, 95% CI [0.450-0.660]) and are referred to as hypomorphic.
Sonic hedgehog signaling (Shh) genes are largely unchanged in kif7 co63/co63 zebrafish: Kif7 is known to function both within sonic hedgehog signaling (Shh) pathway and as an organizer of the axonemal tip of the primary cilia compartment [34]. Common transcriptional targets of hedgehog signaling include ptch1 and gli1 [44,45], which can serve as indicators of Shh pathway activation. We observed that ptch1 transcripts were elevated in both 4 dpf embryos (9.816 fold change, p = 1.16 × 10 −5 ) and 6 wpf scoliotic zebrafish (50.96 fold change, p = 4.02 × 10 −35 ) compared to AB. Gli1 was slightly upregulated in 6 wpf scoliotic kif7 co63/co63 zebrafish compared to AB (1.461 fold change, p =1.05 × 10 −5 ). Additional RT-qPCR data for Shh pathway genes are provided in Supplemental Table S4. Cilia genes as listed in the SYSCILIA database [46] appear to be largely unchanged in kif7 co63/co63 zebrafish as compared to AB (Supplemental Table S4. Although we did not obtain usable data for Shh genes in scoliotic vs. non-scoliotic kif7 co63/co63 zebrafish through RNA sequencing, RT-qPCR results suggested that there was little difference in these genes ( Figure 4, Supplemental Pankeratin staining confirms an upregulation of keratins in scoliotic kif7 co63/co63 zebrafish: To confirm an enrichment of keratins in scoliotic kif7 co63/co63 zebrafish compared to genotype-matched kif7 co63/co63 non-scoliotic zebrafish, we performed pankeratin staining of 6 wpf kif7 co63/co63 zebrafish. Scoliotic zebrafish displayed increased levels of keratins and prekeratins across muscle tissues in whole mount sections ( Figure 5). Additionally, scoliotic zebrafish displayed increased keratins within the intervertebral disc (IVD).

Discussion
In this work, we build on previous findings in which human data supported by a unique zebrafish animal model suggest that mutations in KIF7 may contribute to the IS phenotype [33]. Exploration of our unique zebrafish model, kif7 co63/co63 , which develops spi- Increased keratin/prekeratin can be seen in the IVD of scoliotic fish compared to non-scoliotic fish and wild type, as denoted by the black arrow. N = 3 zebrafish were analyzed for each group (wild type, scoliotic kif7 co63/co63 , and non-scoliotic kif7 co63/co63 ).

Discussion
In this work, we build on previous findings in which human data supported by a unique zebrafish animal model suggest that mutations in KIF7 may contribute to the IS phenotype [33]. Exploration of our unique zebrafish model, kif7 co63/co63 , which develops spinal curvatures during the juvenile period, supports the role of distinct molecular mechanisms which underlie the pathology of the scoliosis phenotype. Our results continue to support the role of hypomorphic mutations in KIF7 as contributory to IS pathogenesis.
As expected, kif7 mRNA was reduced in both scoliotic and non-scoliotic kif7 co63 zebrafish. Overall, zebrafish transcriptomes clustered primarily by sample age and genotype. Regardless of scoliotic phenotype, homozygous kif7 co63/co63 zebrafish clustered together on the PCA plot (Figure 1), indicating that the transcriptomes of these zebrafish were overall similar.
Kif7 as a ciliary kinesin has dual roles within the cell; one, in the sonic hedgehog (Shh) signaling pathway and, second, as a ciliary regulator [34]. Kif7 is broadly expressed in human [47] and zebrafish [48]. In mouse, Kif7 functions as both a negative and positive regulator of the Shh pathway within the primary cilium downstream of Smoothened (Smo) and upstream of Gli transcription factors [49]. Ptch1, which encodes Patched, the primary receptor for Shh, is noted to be upregulated in 4 day (9.816 fold change, p = 1.16 × 10 −05 ) and 6 wpf scoliotic kif7 co63/co63 zebrafish (50.96 fold change, p = 4.02 × 10 −35 ) as compared to age-matched AB zebrafish. Although we did not obtain usable bulk RNAseq data for Shh genes in scoliotic vs. non-scoliotic kif7 co63/co63 zebrafish, RT-qPCR data showed little difference in Shh genes in 6 wpf scoliotic vs. non-scoliotic kif7 co63/co63 zebrafish, with the exception of dlg5a (15.225-fold change, 95% CI [9.727, 23.830]). Dlg5a binds directly to Kif7 [50], and was observed as upregulated in kif7 co63/co63 embryos over controls within our previous study [33]. Dlg5 is also a regulator of Gli1 protein ubiquitination and degradation and, most recently, is recognized as a major factor of Shh signaling, particularly in relation to glioblastoma tumors, reflecting its known role in critical processes involving cell-cell adhesion during neural development [51,52]. Although we cannot rule out the possibility that Shh signaling plays a tissue-specific role in spine morphogenesis in kif7 co63/co63 zebrafish, as Shh target genes ptch1 and gli1-2 were largely unchanged between scoliotic and non-scoliotic mutants, the upregulation of dlg5 could represent alternative functional mechanisms of cellular interactions.
We observed striking differences in several other transcripts between scoliotic and nonscoliotic kif7 co63/co63 zebrafish, notably a large upregulation in the zebrafish cytoskeletal keratins krt4 and zgc:158846 (krtt2c6), which are orthologous to human KRT8 (fold change= 174.85, p = 2.79 × 10 −289 ) and KRT7, −8 and −9 (fold change = 107.63, p = 1.37 × 10 −159 ). This upregulation of zebrafish keratin krt4 was too variable between specimens to draw a definitive conclusion by RT-qPCR, which may have been due to a difference in the methodologies used [53]. However, we also observed increased levels of keratins and prekeratin proteins in the dermis, musculature, and intervertebral disc (IVD) of scoliotic kif7 co63/co63 zebrafish as compared to non-scoliotic genotype-and agematched mutant zebrafish. These key transcriptomic differences may arise spontaneously and contribute to the development of spinal curvatures in scoliotic mutants, or they may be a secondary effect that develops after the initiation of spinal curvatures.
It is unknown specifically how kif7 mutations may lead to increased expression of keratins, a known product of keratinocytes, in hypomorphic kif7 co63/co63 zebrafish. Due to the role of Shh signaling in the pathogenesis of basal cell carcinoma of the skin, the role of Kif7 has been studied in keratinocytes. This work has been primarily performed in the mouse as the model system; however, studies within the zebrafish have shown similarities of epidermal ontogeny [54,55]. In keratinocytes, Kif7 has been shown to have a crucial role as a positive regulator of Shh signaling through its regulatory function of the Gli transcription factors [56]. In embryonic keratinocytes, inactivation of both Kif7 and Sufu, an additional regulator of Gli transcription factors, leads to a loss of epidermal differentiation and follicular fate, while in the adult this loss leads to the induction of basal cell carcinoma [57]. As seen in histological sectioning (Figure 4), our hypomorphic kif7 co63/co63 zebrafish with scoliosis exhibited significant upregulation of keratin expression overall. The identification of the specific mechanism as to what drives this strong expression of keratin within the muscle and epithelium of the hypomorphic scoliotic kif7 co63/co63 zebrafish remains to be determined.
As stated, the specific human orthologues of the top differentially expressed genes within our expression data are KRT7 and KRT8. KRT7 and KRT8 encode keratin-7 and keratin-8, respectively, and are members of the type II keratin family that are broadly expressed and form intermediate filaments in the cytoplasm of epithelial cells. KRT7 expression has a role in epithelial-mesenchymal transition and its dysregulation has been highly associated with various types of cancers and tumor progression [58][59][60]. Cytokeratins of the intermediate filament subgroup, which includes KRT8, have been observed in adult mouse skeletal muscle [61]. Keratin-8, in addition to keratin-19, links the contractile apparatus of striated muscle to dystrophin [62]. More recently, KRT8, KRT18 and KRT19, all recognized markers of notochordal cells, have been found to be expressed in a subpopulation of adult nucleus pulposus (NP) cells. NP cell degeneration is associated with intervertebral disc degeneration (IVDD), a prominent health problem worldwide [63,64]. The finding of these cellular markers within the NP has allowed for study of the intervertebral disc (IVD) and its degeneration in both zebrafish and humans [65]. The role of keratins within the IVD as a potential molecular mechanism associated with the onset of scoliosis warrants further study. Future explorations will focus on isolating the tissue or cell type driving this expression pattern.
Ces2 and TTC38 were significantly downregulated within the scoliotic homozygous kif7 co63/co63 zebrafish when compared to non-scoliotic homozygous kif7 co63/co63 fish. Both genes are related to metabolic abnormalities through enzymatic reactions of various drugs and endogenous compounds (via NCBI RefSeq [66]), to which there is no known relationship to axial skeletal development.
Collectively, we present transcriptomic data that associate an upregulation of cytokeratins with the phenotype of scoliosis within a kif7 co63/co63 zebrafish animal model. This comparative transcriptomics analysis has significance in the determination of unique differentially expressed genes as related to the scoliosis phenotype. This initial step, when followed with gene-set-enrichment analyses and orthology mapping, has the potential to reveal significant associations with biological relevance to the observed phenotype of idiopathic scoliosis.